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Abstract 

We define Wannier functions for interacting systems, and show that the results on 
the localization of the Wannier functions for non-interacting systems carry over to 
the Wannier functions for interacting systems. In addition we demonstrate that the 
characterization of metals and insulators by the decay properties of their respective 
density matrices does not only apply to non-interacting, but also to interacting 
systems. As a prototypical example of a correlated system we investigate the one- 
dimensional Hubbard model. We propose an expression for the density matrix of 
that model, and derive a relation between the decay constant of the density matrix 
and the gap. 
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The locality properties of solids, their so called nearsightedness, have recently 
moved to the focus of much attention [1-7]. This is largely due to the intense 
efforts going into the development of electronic structure methods that scale 
linearly with system size [5,6]. These methods depend crucially on the locality 
of the density matrix. The closely related issue of the locality of the Wannier 
functions [8] has also attracted renewed interest, largely due to the develop- 
ment of methods for the practical construction of localized Wannier functions 
[4]. The uses of Wannier functions are well proven. First, they are fundamental 
in the theory of electron dynamics in the presence of weak external fields [9] . 
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Second, they allow for an intuitive interpretation of the bonding properties in 
solids [4]: Localized Wannicr functions correspond to cither bonds or lone elec- 
tron pairs. Third, they arc at the center of the modern theory of polarization 
[10]. Finally, they are important for some linear scaling algorithms [5,6]. 

General results on the localization of the Wannier functions are quite difficult 
to obtain. The problem has actually been called "one of the few basic questions 
of the quantum theory of periodic solids in the one-electron approximation 
which is not completely solved" [11]. The most important results so far are 
the proofs of the existence of exponentially localized Wannier functions for 
isolated, simple bands in any dimension [11], and for complex bands in the 
tight-binding limit, and in perturbation theory [12]. For practical applications 
knowledge of the decay constants is of considerable importance. First results 
were obtained in [13], considerably more general results for the density matrix 
were given in [2]. 

Trying to extend the results to interacting systems we are faced with a fun- 
damental problem: The Hamiltonian is no longer a one-body operator. Hence 
it seems to be impossible to define meaningful single-electron Bloch waves 
and the corresponding Wannier functions. This is true already in the limit of 
weak interactions, although there one would expect that the non-interacting 
Wannier functions still could be useful. We can, however, avoid the problem of 
having to deal with a many-body operator, by changing perspective. Looking 
at the one-body density matrix, we realize that it can replace the (single- 
electron) Hamiltonian in the standard construction of Wannier functions. We 
are thus led to Wannicr functions made from natural orbitals [15], which we 
christen natural Wannicr functions. Instead of the energy bands ej(k) we now 
consider occupation bands nj{]i). Wc prove that the projection operator P(k) 
onto an isolated set of occupation bands is analytic in k. Once this result is 
estabhshed, the results [11,12] obtained for the standard Wannier functions 
immediately carry over to the natural Wannier functions. So the natural Wan- 
nier functions can be considered as the natural generalization of the concept 
of Wannier functions to interacting systems. Being constructed from natural 
orbitals, they are an optimal basis, meaning that it is sufficient to consider 
only the bands nj(k) with high occupation to obtain a good description of 
the interacting system [15]. Being Wannier functions, they are localized, thus 
allowing to take advantage of 0{N) methods in evaluating, e.g.. Coulomb 
matrix elements. Moreover, they can be expected to represent the chemical 
bonding in the correlated system, as do the ordinary Wannier functions in the 
independent-particle case. We furthermore show that the characterization of 
metals and insulators in terms of the decay properties of the density-matrix 
[14], also applies to correlated systems. Thus the decay constant 7 may play 
a similar role as the localization length defined in [7]. Finally, we consider the 
one-dimensional Hubbard model as an explicit example of a correlated system, 
deriving a relation giving the decay constant 7 as a function of the gap, thus 
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extending the results of [2] to an interacting system. 

We start our construction from tlie many-body wavefunction \l/(ri, r2, Ttv). 
For simplicity we consider spinless electrons. The one-body density matrix D 
is then given by 

^(r, r') = [ dr2--- (ir^^'*(r', rs, . . . r^^)^'(r, rs, . . . r^^). 



Its eigenfunctions are the natural orbitals, its eigenvalues the natural occupa- 
tion numbers. Clearly, D is Hermitian. Furthermore, translating all the spatial 
arguments of ^ by a lattice vector R multiplies the wavefunction by a phase 
factor. Since these phase factors cancel inside the density matrix, we find 
D(r + R, r' + R) = D{r, r'). Thus the Bloch theorem applies and, in analogy 
to the Hamiltonian of a periodic solid in an independent-particle picture [16], 
the eigenvalues of the density matrix form bands rij (k) , the occupation bands, 
with the corresponding natural orbitals being Bloch functions 

$,,k(r) = e^^^ C/,,k(r) , (1) 



where j is the band index, k a vector in the Brillouin zone, and the f/^ ^(r) 
are periodic functions with respect to the real space primitive cell. They are 
the eigenvectors of the k-dependent density matrix 

Dk = E^.(k)|t^.,k)(t/,-,k|. 
j 



Using the standard prescription [16], we can then construct Wannier functions, 
which, for obvious reasons, we call natural Wannier functions: 

^ ^ BZ 



here V is the volume of the real-space cell and the integration is over the 
Brillouin zone (BZ). 

The construction of the natural Wannier functions is very similar to that of 
the conventional Wannier functions for a non-interacting system, the density 
matrix D taking the place of the non-interacting Hamiltonian. In the limit 
of vanishing interaction, the natural Wannier functions do, however, not re- 
duce to the conventional Wannier functions. This is clear, since in the limit 
of no interaction the density matrix becomes the projector onto the occupied 
subspace; i.e. all occupied (unoccupied) states are degenerate with eigenvalue 
one (zero). By degenerate perturbation theory, in the limit of vanishing inter- 
action, the natural Wannier functions will therefore have to diagonalize the 
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first term in the perturbation expansion in their respective subspace (occupied 
or unoccupied) [15]. They are therefore gencrahzcd Wannier functions in the 
sense of [4], which, respectively, span the space of occupied and unoccupied 
states. Clearly, using several bands increases the fiexibihty to construct more 
localized Wannier functions, since one can take advantage of the unitary trans- 
formations allowed in the space spanned by the Bloch functions of different 
bands. In the extreme limit, where we allow all bands, it is obviously possible 
to construct Wannier functions that are perfectly localized delta functions. In 
this respect it is interesting to note that quantum Monte Carlo calculations 
of the lowest natural orbitals show that they are very similar to the occupied 
Kohn-Sham orbitals [17]. Hence the corresponding natural Wannier functions 
are also expected to be very similar to their Kohn-Sham counterparts. This 
similarity is presumably no longer valid for the Wannier functions arising from 
the unoccupied bands. Whereas the virtual Kohn-Sham Wannier functions are 
much less localized than the occupied ones [18] this seems not to be the case 
for the natural Wannier functions. 

We now want to show under what conditions the natural Wannier functions 
are exponentially localized. The key input to the corresponding proofs for 
the ordinary Wannier functions [11,12] is the analyticity of the projection 
operator onto the band states. The proof for energy-bands was given in [14]. 
For the natural Wannier functions we have to prove the analyticity on a strip 
K = {k' -|- ih", |k"| < A} of the projector P(k) onto an isolated set B of 
occupation bands. Because of the special properties of the density matrix, the 
proof is straightforward. is hermitian for real k; A^-representability [15] 
requires for the eigenvalues rii G [0,1]; moreover, since a unit cell contains 
only a finite number of electrons and the bands are continuous, TrDk is finite, 
i.e. is trace class, and therefore, in particular, compact. If we assume that 
Dk is analytic on K, we can apply the analytic Predholm theorem [19], which 
guarantees that the resolvent of Dk is meromorphic on K with the residues at 
the eigenvalues being finite rank operators. Choosing a contour C in K, which 
exclusively encircles all occupation numbers in B (here we use that the bands 
in B are isolated, i.e. they do not intersect with, or touch any band not in B), 
we obtain the occupation band projector 



which, as desired, is analytic on K. Given the analytic band projection opera- 
tor, the proofs given in [11,12], with the Hamiltonian replaced by the one-body 

density matrix, guarantee the existence of exponentially localized Wannier 
functions if B contains only a single band, or, for complex bands B, in the 
tight-binding limit and in perturbation theory (of Dk) around a situation, 
where there do exist exponentially localized Wannier functions. 




(3) 
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In the proof given above, we assumed that the k-dependent density matrix 
is analytic in k. It is then natural to ask what the analyticity of means. 
Let us therefore distinguish two cases: Either is analytic in k or it has 
a non-analyticity at some k = kp. Then, in the absence of degeneracies, in 
the first case, the occupation band structure nj(k) will be analytic in k [19], 
while in the second case it will have a non-analyticity at kp. In the first case 
there is no Fermi surface, and we thus associate this case with an insulator. 
As for the second case, assuming that perturbation theory holds, a discon- 
tinuity in nj{k) implies a Fermi liquid, while an algebraic singularity would 
point to Luttinger liquid behavior. We thus associate this case with a metal. 
But we note that we are not aware of a proof that a discontinuity in nj (k) 
by itself guarantees metallic behavior even when perturbation theory breaks 
down [20]. It is, however, hard to imagine an insulator with such a discontinu- 
ity, in particular, since it would also show up in the momentum distribution 
-^(p) — I drdr'exp(— ip(r — r'))£)(r, r'), which is, in contrast to ^^(k), an 
experimentally accessible quantity. From the Paley- Wiener theorem [19] (see 
also [14]), it follows that for an insulator the density matrix, being the Fourier 
transform of Dy_{r,r'), decays exponentially with increasing distance |r — r'|. 
In the three-dimensional case, the decay constant will in general be different 
along different directions. For a metal, on the other hand, where £)k(i"',r) is 
non-analytic on the Fermi surface, the density matrix decays only algebraically. 
The decay properties of the density matrix for interacting systems are thus 
qualitatively the same as for non-interacting systems [14]. We emphasize that 
this result is valid for zero temperature, while at finite T the decay should 
become exponential also for a metal [21]. 

As a specific and important example, we now analyze the localization prop- 
erties of the density- matrix for a prototypical correlated system, the one- 
dimensional Hubbard model (with lattice constant a and nearest-neighbor 
hopping matrix element t) [22] 

H^-t Yl (iaCia + UY1 ^n^a ■ (4) 

For U — Q the system is metaUic, while for any finite C/ it is a Mott insulator 
[23]. Motivated by the result shown in Fig. 1 we choose as an Ansatz the 
product of the exact density matrix for [/ = and an exponential factor: 

D„,„=(*|cLc„,|*) = 2i!^<^e--H. (5) 



To estimate the decay constant, we calculate Dq i from the exact ground-state 
energy of the Hubbard chain [23]. Using the Hellmann-Feynman theorem we 
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Fig. 1. Exponential decay of the density matrix for the Hubbard chain. The symbols 
give the numerical values of the density matrix for a chain of 16 sites. The lines 
connect the values of Ansatz (5), using the decay constants given in (7). 

find 



(6) 



which, together with the Ansatz (5), yields 



oo 



T<. = -ln|.yd. ^°W^4P |. (7) 



where Jq and Ji arc Bessel functions. Fig. 2 shows a plot of 7a as a function 
of the gap Eg/t = U/t-A + S j'^ dxJi{x)/{x{l + ex^{Ux/2t))). For large Eg 
we find the asymptotic behavior 



while for small gap there exists no expansion, as for the Hubbard chain the 
point = is non-analytic. The behavior of 7a is thus qualitatively different 
in both the large and small gap limit from the analytical results for non- 
interacting systems (cf. Fig. 1 in Ref. [2]), although the overall shapes of 
the curves look similar. We note that the decay of the density matrix (5) 
has a power-law exponent different from the universal exponent for the non- 
interacting density- matrix found in [24]. 

Fourier transforming the Ansatz (5), we obtain 

,,,11 / cos(A;a) \ , . 

nik) = - + - arctan — tt—^ ■ (9) 
2 TT \smh(7a) J 
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Fig. 2. Decay constant 7a of the density matrix for the Hubbard chain as function 
of the gap Eg. The dashed hne shows the large-?/ approximation (8). 
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Fig. 3. Occupation band structure for periodic Hubbard chains (determined by 
exact diagonahzation) compared to (9) with 7a obtained from (7). Ah curves stih 
differ significantly from the expression for n{k) in the large-?/ limit [25], which, for 
U = 12, is plotted for comparison. 

By construction, this is the step function for [/ — > 0, while in the hmit of 
large U wc have n{k) = | + ^^^^j^ cos(A;a); i.e. our Ansatz (5) is exact, both 
for U ^ and U ^ 00 [25]. For intermediate U, we compare (9) with the 
momentum distribution obtained from exact diagonahzation of finite Hubbard 
chains. This is shown in Fig. 3: For U large, but still far from the large-C/ limit, 
the agreement is perfect, and even for fairly small U, where the decay length 
1/70 is of the order of the chain length, the agreement between (9), which was 
derived for the infinite chain, and the results for the finite rings is amazingly 
good. 

The half- filled Hubbard chain is a Mott insulator for [/ > 0, only for [/ = 
it is a metal. As expected, this is reflected in the decay properties of the 
density matrix. While for f/ = it decays as 1/m, for any finite U the decay 
is exponential. It is interesting to compare the decay constant 7a with the 
localization length A/a = \/d/27r, with d — — Huin^oo N In \zn\'^ and zn — 



7 



(^ole^^'^'l^o), as defined in [7]. We find that both criteria describe the 
metal-insulator transition correctly. As a surprising fact we note that 1/A 
seems to be linear in U down to very small values of U, with the constant of 
proportionality given by a large U expansion. 

Away from half-filling we have tried an Ansatz in the spirit of (5) of the form 
smikpam) |^ _|_ j^j-o!^ Fouricr transforming leads to an n(k) with a Luttinger- 
like singularity with exponent a at kp- It fails, however, to also produce a 
singularity at 3k p [26]. 

To summarize, we have defined natural Wannier functions for interacting sys- 
tems from the eigenfunctions of the density matrix, and have shown under 
what conditions they can be proven to be exponentially localized. The natural 
Wannier functions provide an optimal, localized basis for describing a corre- 
lated system. In terms of computational efficiency, they will allow the use of 
0{n) methods in many-body calculations. We also expect them to provide 
understanding of the bonding in correlated solids. In addition, we have shown 
that the characterization of metals and insulators by the decay properties of 
the density matrix does also apply to interacting systems, at least as long as 
perturbation theory holds. Finally, we have investigated the one-dimensional 
Hubbard model, proposed an expression for the density matrix of that model, 
and, for this Mott insulator, derived a relation between the decay constant of 
the density matrix and the gap. 

We thank L.N. Trefethen for helping in identifying the relevant mathematical 
theorems, and T. Arias, O. Gunnarsson, P. Horsch, J. Hutter, R.M. Martin, 
G. StoUhoff, and D. Vanderbilt for interesting discussions. 

Note added in proof: After finishing the present work we became aware of 
the paper [27], in which also introduces natural Wannier functions. We note, 
however, that in that paper the analyticity of the Bloch functions for simple 
bands is take for granted, rather than proven. 
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